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Sea-level  rise  during  the  last  deglaciation  and  through  the  Holocene  was  influenced  by  deformational, 
gravitational,  and  rotational  effects  (henceforth  glacial  isostatic  adjustment,  GIA)  that  led  to  regional 
departures  from  eustasy.  Deglacial  sea-level  rise  was  particularly  variable  spatially  in  areas  adjacent  to 
the  Cordilleran  and  Laurentide  Ice  Sheets.  Such  regional  variability  in  sea  level  due  to  GIA  is  important  to 
identify  when  investigating  potential  coastal  migration  pathways  used  by  early  Americans.  An  improved 
understanding  of  regional  sea-level  rise  may  also  be  used  for  predictive  modeling  of  potential  archae¬ 
ological  sites  that  are  now  submerged.  Here  we  compute  relative  sea-level  change  across  the  California 
—Oregon— Washington  and  Bering  Sea  continental  shelves  since  the  Last  Glacial  Maximum  using  an  ice- 
age  sea-level  theory  that  accurately  incorporates  time-varying  shoreline  geometry.  The  corresponding 
non-uniform  sea-level  rise  across  these  continental  shelves  reveals  significant  departures  from  eustasy, 
which  has  important  implications  for  improved  understanding  of  potential  coastal  migration  routes  and 
predictive  modeling  of  the  location  of  now-submerged  archaeological  sites. 

Published  by  Elsevier  Ltd. 


1.  Introduction 

Following  the  Last  Glacial  Maximum  (-26-19  ka)  (ka:  kilo- 
annum.  All  years  reported  herein  are  in  calendar  years)  (Clark  et  al., 
2009),  ice-sheet  retreat  caused  globally  averaged  sea  level  to  rise  by 
-130  m  (Austermann  et  alM  2013;  Yokoyama  et  al.,  2000),  with 
present  sea  level  reached  -5  ka  (Woodroffe  et  al.,  2012).  During  this 
period  of  lower  sea  level,  large  expanses  of  continental  shelves 
were  exposed,  providing  coastal  migration  routes  and  occupation 
sites  that  have  since  become  submerged  as  global  sea  level  rose  to 
its  present  height.  Reconstructing  the  paleogeography  of  these 
now-submerged  landscapes  is  thus  important  for  inferring  the 
most  likely  routes  taken  by  early  coastal  migrations,  as  well  as  for 
predicting  the  locations  of  subsequent  occupation  sites. 

These  issues  have  recently  become  of  wide  interest  along  the 
west  coast  of  the  Americas,  where  new  evidence  for  pre-Clovis 
cultures  and  early  maritime  activity  has  focused  attention  on  the 
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importance  of  coastal  paleogeography,  particularly  with  regard  to 
the  peopling  of  the  Americas  (Graf  et  al.,  2013).  For  much  of  the 
20th  century,  the  prevailing  paradigm,  referred  to  as  the  “Clovis 
First”  model,  posited  that  the  initial  colonization  of  the  Americas 
south  of  the  Laurentide  and  Cordilleran  ice  sheets  was  by  migration 
from  Beringia  through  an  ice-free  corridor  that  developed  between 
the  two  ice  sheets  during  the  last  deglaciation.  The  model  was  first 
developed  because  the  locations  of  the  first  Clovis  sites  were  south 
of  the  ice-sheet  corridor,  suggesting  an  initial  entry  point.  The 
paradigm  was  subsequently  supported  by  the  evidence  that  the  age 
for  the  opening  of  the  ice-free  corridor  (Catto,  1996)  roughly 
coincided  with  the  earliest  age  of  Clovis  artifacts  (Haynes,  2002), 
with  the  attendant  implication  that  prior  closure  of  the  corridor 
would  have  prevented  earlier  migrations.  The  model  implied  that 
upon  arriving  to  the  interior  of  the  continent  south  of  the  ice  sheets, 
these  early  peoples  slowly  migrated  towards  coastal  regions,  where 
they  adapted  to  the  local  environment  (Erlandson  et  al.,  2008). 

Despite  the  prevailing  view  of  the  Clovis  First  model,  several 
questions  were  raised  as  to  whether  humans  would  be  able  to 
survive  the  harsh  environmental  conditions  and  biologically  un¬ 
productive  landscape  of  the  ice-free  corridor  (Fladmark,  1979; 
Mandryk,  1993;  Meltzer,  2004).  Largely  for  these  reasons, 
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Fladmark  (1979),  building  on  an  earlier  suggestion  by  Heusser 
(1960),  proposed  an  alternative  migration  route  along  the  west¬ 
ern  coast  of  North  America,  pointing  out  that  a  “chain  of  refugia” 
along  the  coast  would  have  been  more  environmentally  suitable  for 
occupation,  given  the  available  marine  resources  and  the  possibility 
that  the  early  people  had  simple  watercraft  for  navigating  the  route. 

A  number  of  developments  in  the  last  ten  years  have  further 
challenged  the  Clovis  First  model  to  the  extent  that  some  have  now 
suggested  that  it  is  all  but  defunct  (Adovasio  and  Pedler,  2013; 
Erlandson,  2013).  The  demise  of  the  model  has  been  largely 
driven  by  mitochondrial  DNA  (Fagundes  et  al.,  2008;  Perego  et  al., 
2009)  and  archeological  (Dillehay  et  al.,  2008;  Gilbert  et  al.,  2008; 
Goebel  et  al.,  2008;  Jenkins  et  al.,  2012;  Waters  et  al.,  2011a; 
Waters  and  Stafford,  2007;  Waters  et  al.,  2011b)  studies  which 
provide  compelling  evidence  for  pre-Clovis  occupation  of  the 
Americas  by  14-15  ka. 

Based  on  current  understanding,  this  earlier  age  for  the  arrival  of 
First  Americans  south  of  the  ice  sheets  likely  preceded  the  opening 
of  the  ice-free  corridor  by  several  hundred  to  a  thousand  years 
(Dyke,  2004),  although  dating  of  post-glacial  eolian  sediments  with 
optically  stimulated  luminescence  suggests  the  possibility  that  the 
corridor  was  open  early  enough  for  a  migration  route  by  pre-Clovis 
peoples  (Munyikwa  et  al.,  2011).  If  the  ice-free  corridor  route  is  no 
longer  viable,  however,  then  by  default  a  coastal  route  is  required, 
with  a  route  along  the  west  coast  of  North  America  now  considered 
the  most  likely  (Erlandson,  2013).  Several  lines  of  evidence  have 
bolstered  arguments  for  this  route,  including  additional  arguments 
for  the  wealth  of  resources  that  the  route  offered  (Erlandson  et  al., 
2007),  the  older  ages  of  human  settlement  and  maritime  activity 
along  the  Pacific  Coast  (Erlandson  and  Braje,  2011 ;  Erlandson  et  al., 
2008,  2011),  and  the  identification  of  areas  along  the  Alaskan  and 
British  Columbian  coasts  that  were  habitable  one  to  three  millennia 
before  the  opening  of  the  ice-free  corridor  (Mandryk  et  al.,  2001; 
Mann  and  Flamilton,  1995;  Misarti  et  al.,  2012). 

Reconstructing  the  paleogeography  that  existed  at  the  time  of  a 
coastal  migration,  as  well  as  during  subsequent  times  when  peo¬ 
ple  inhabited  now-submerged  shelf  areas,  requires  knowing  the 
sea-level  history  during  the  last  20,000  years.  A  common  approach 
to  reconstruct  coastal  paleogeography  (Anderson  et  al.,  2013, 
2010;  Davis  et  al.,  2004;  Kennett  et  al.,  2008;  Manley,  2002)  is 
to  uniformly  lower  sea  level  by  the  amount  suggested  from  far- 
field  (e.g.,  Barbados,  Tahiti)  sea-level  records  (Bard  et  al.,  1996; 
Fairbanks,  1989)  which  are  assumed  to  represent  the  global 
mean.  This  uniform,  global  mean  sea-level  change  is  commonly 
referred  to  as  the  “eustatic”  change.  However,  adopting  the 
eustatic  change  to  reconstruct  coastal  paleogeography  moves 
beyond  this  simple  definition  to  an  implicit  assumption  that 
glacial  meltwater  entered  the  oceans  in  a  geographically  uniform 
manner  (i.e.,  the  so-called  bathtub  model).  Sea-level  rise  during 
the  last  deglaciation,  however,  was  influenced  by  glacial  isostatic 
adjustment  (GIA)  associated  with  the  exchange  of  mass  between 
ice  sheets  and  oceans  that,  for  the  majority  of  the  ocean  basins,  led 
to  significant  regional  departures  from  eustasy  (Milne  and 
Mitrovica,  2008). 

Deglacial  sea-level  rise  across  the  California- 
Oregon-Washington  and  Bering  Sea  continental  shelves  would 
have  been  particularly  variable  spatially  as  a  consequence  of  their 
relative  proximity  to  the  Cordilleran  and  Laurentide  Ice  Sheets.  In 
this  article,  we  use  a  state-of-the-art  numerical  model  of  post¬ 
glacial  sea-level  change  (Kendall  et  al.,  2005)  to  predict  relative 
sea-level  (RSL)  change  in  the  region  of  the  Channel  Islands,  Cal¬ 
ifornia,  and  across  the  Oregon-Washington  and  Bering  Sea  con¬ 
tinental  shelves.  We  compare  these  predictions  to  sea-level 
changes  that  would  be  associated  with  a  globally  uniform 
“eustatic”  rise. 


2.  Sea-level  modeling 

The  ice-age  sea-level  calculations  were  performed  using  an  al¬ 
gorithm  (Kendall  et  al.,  2005;  Mitrovica  and  Milne,  2003)  that  re¬ 
quires,  as  input,  models  for  both  the  global  geometry  of  ice  sheets 
over  the  last  glacial  cycle  and  the  Earth's  viscoelastic  structure  (see 
below).  The  algorithm  yields  a  gravitationally  self-consistent  ocean 
redistribution  and  it  accounts  for  the  viscoelastic  deformation  of 
the  solid  Earth  in  response  to  the  (ice  plus  water)  loading,  as  well  as 
associated  perturbations  to  the  Earth's  gravitational  field  and 
rotational  state  (Mitrovica  et  al.,  2005).  The  sea-level  theory  accu¬ 
rately  treats  for  the  migration  of  shorelines  due  to  local  sea-level 
variations  and  changes  in  the  perimeter  of  grounded,  marine- 
based  ice  sheets  (Milne  et  al.,  2002),  and  an  iterative  procedure 
adopted  within  the  algorithm  guarantees  that  the  predicted 
present-day  topography  matches  the  observed  topography 
(Kendall  et  al.,  2005).  Sediment  redistribution,  including  both 
erosion  and  deposition,  are  not  included  in  the  modeling. 

All  our  GIA  calculations  are  based  on  a  pseudo-spectral  solver 
truncated  at  spherical  harmonic  degree  256,  which  provides  a 
surface  spatial  resolution  of  GIA  effects  to  a  level  of  approximately 
100  km  (Kendall  et  al.,  2005).  In  results  presented  below,  these 
computed  changes  in  global  sea  level  are  superimposed  onto  a 
high-resolution  regional  grid  of  modern  topography  to  track 
changes  in  shoreline  position  as  a  function  of  time.  Since  GIA- 
induced  changes  in  sea  level  have  relatively  smooth  spatial  geom¬ 
etries,  adopting  a  higher  truncation  in  the  pseudo-spectral  solver 
would  not  significantly  alter  the  predicted  evolution  of  shoreline 
positions.  In  the  calculations  presented  below,  we  used  the  3- 
arcsecond  U.S.  Costal  Relief  Model  (http;// www.ngdc.noaa.go v/ 
mgg/coastal/crm.html)  for  CA-OR-WA  and  the  1-arcminute 
ETOPOOl  topography  grid  for  Beringia  (CRM  coverage  is  not  avail¬ 
able  for  this  region)  (Amante  and  Eakins,  2009). 

We  adopt  a  spherically  symmetric,  self-gravitating,  linear 
(Maxwell)  viscoelastic  Earth  model.  The  elastic  and  density  struc¬ 
ture  of  this  model  is  prescribed  from  the  seismic  model  PREM 
(Dziewonski  and  Anderson,  1981).  In  an  earlier,  preliminary  study 
(Clark  et  al.,  2014)  of  shoreline  migration  in  the  same  region,  we 
adopted  the  so-called  VM2  radial  profile  of  mantle  viscosity,  which 
is  coupled  to  the  ICE-5G  ice  history  since  the  last  interglacial 
(Peltier,  2004).  The  VM2  model  (Peltier,  2004)  is  characterized  by  a 
relatively  moderate,  factor  of  5,  increase  in  viscosity  from  the  base 
of  a  90  km  thick,  high  viscosity  (effectively  elastic)  lithosphere  to 
the  core-mantle-boundary.  In  contrast,  the  simulations  described 
below  adopt  a  model  of  ice  history  developed  at  the  Australian 
National  University  (Fleming  and  Lambeck,  2004)  which  is  coupled 
to  an  Earth  model  with  a  lithospheric  thickness  of  96  km,  an  upper 
mantle  viscosity  of  0.5  x  1021  Pa  s,  and  a  lower  mantle  viscosity  of 
8  x  1021  Pa  s.  This  Earth  model  (henceforth  the  LM  model)  is 
consistent  with  inferences  from  a  number  of  studies  based  on 
globally  distributed  sea-level  data  sets  (Mitrovica  and  Forte,  2004; 
Nakada  and  Lambeck,  1989).  Our  adoption  of  the  LM  model  is  also 
motivated  by  the  study  of  Muhs  et  al.  (2012),  who  examined  the 
sea-level  history  over  the  last  glacial  cycle  at  San  Nicolas  Island,  CA, 
one  of  the  Channel  Islands,  with  emphasis  on  sea-level  oscillations 
through  marine  isotope  stage  (MIS)  5.  They  found  that  MIS5a  and 
5c  highstands  at  San  Nicolas  Island,  as  well  as  previously  published 
highstands  at  Barbados  and  the  Florida  Keys,  were  best  fit  by  a  GIA 
prediction  based  on  a  model  similar  to  our  LM  model,  but  were  not 
well  fit  by  the  model  VM2. 

In  results  described  below  we  will  also  consider  the  sensitivity 
of  some  of  our  predictions  to  the  presence  of  lateral  variations  in 
Earth  structure,  including  mantle  viscosity  and  lithospheric  thick¬ 
ness.  We  also  incorporate  tectonic  plate  boundaries  as  zones  of 
particularly  low  viscosity.  These  calculations  are  based  on  a  finite 
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Age  (ka) 

Fig.  1.  Predicted  relative  sea  level  (RSL)  curves  (blue)  for  two  sites  within  the  Channel 
Islands  region  compared  to  the  eustatic  curve  for  the  ANU  deglaciation  history  (red). 
The  RSL  curves  are  computed  for  sites:  34°N  120°W  (solid  blue  line)  and  33.5°N  118°W 
(dashed  blue  line).  (For  interpretation  of  the  references  to  color  in  this  figure  legend, 
the  reader  is  referred  to  the  web  version  of  this  article.) 

volume  numerical  model  of  ice  age  adjustments  (Latychev  et  al., 
2005). 

3.  Channel  Islands,  California 

Maritime  peoples  occupied  California's  Channel  Islands  by  at 
least  ~13  ka,  with  the  implication  that  they  used  seaworthy  boats  to 
travel  to  the  islands  from  the  mainland  (Erlandson  et  al.,  2005, 
1996,  2011;  Rick  et  al.,  2001,  2005).  Kennett  et  al.  (2008)  adopted 
a  eustatic  history  inferred  from  far-held  sea-level  records  (Bard 
et  al.,  1996;  Fairbanks,  1989)  to  reconstruct  the  islands'  changing 
paleogeographies  and  distance  to  the  main  coastline  associated 
with  postglacial  sea-level  rise,  with  attendant  implications  for 
changes  in  local  marine  resources.  Reeder-Myers  et  al.  (submitted 
for  publication)  are  also  considering  GIA-induced  changes  in 
shoreline  position  in  the  vicinity  of  Channel  Islands  with  a  focus  on 
ecological  changes  and  human  history. 


We  first  show  examples  of  two  RSL  curves  from  sites  within  the 
Channel  Islands  region,  and  compare  these  to  the  eustatic  curve  for 
the  ANU  ice  history  adopted  in  the  numerical  simulation  (Fig.  1). 
The  RSL  curves  show  a  marked  departure  from  an  assumption  of 
eustasy.  Near  the  end  of  the  LGM  at  ~21  ka,  RSL  was  -25-30  m 
higher  than  would  have  been  predicted  by  assuming  eustatic  sea 
level.  For  periods  following  -13  ka,  the  sign  of  the  discrepancy 
changes  and  the  predicted  RSL  is  as  much  as  -15  m  lower  than 
eustatic  curve.  Note  that  the  multi-meter-scale  differences  between 
the  predicted  and  eustatic  curves  continue  up  to  the  present,  which 
has  implications  for  predicting  sites  in  present-day  estuarine  re¬ 
gions  that  would  have  remained  exposed  for  much  longer  than 
otherwise  predicted  from  the  eustatic  assumption. 

The  Channel  Islands  cover  a  large  enough  area  that  the  predicted 
GIA  contribution  varies  non-negligibly  across  the  area.  As  an 
example,  at  the  LGM,  the  predicted  RSL  at  Santa  Rosa  Island 
(120°W,  34°N)  is  -5  m  lower  than  the  prediction  at  the  present  Los 
Angeles  coastline  (118°W,  33.5°N)  (Fig.  1).  This  difference  reflects 
spatial  gradients  in  the  predicted  RSL  that  are  highlighted  in  Fig.  2, 
which  shows  regional  maps  of  RSL  from  20  ka  to  6  ka  predicted 
using  the  gravitationally  self-consistent  sea-level  theory.  The  pre¬ 
dicted  RSL  on  these  maps  is  more  negative  toward  the  southwest 
(i.e.,  sites  toward  the  southwest  would  have  experienced  greater 
post-LGM  sea-level  rise),  reflecting  the  combined  effects  of  crustal 
subsidence  within  the  peripheral  bulge  of  the  Laurentide- 
Cordilleran  Ice  Sheets  and  migration  of  water  away  from  these  ice 
sheets  as  they  lost  mass  and  thus  exerted  less  of  a  gravitational  pull 
on  the  oceans  (Milne  and  Mitrovica,  2008).  The  predicted  gradient 
is  largest  during  the  LGM  and  diminishes  with  time  (Figs.  1  and  2). 

Fig.  3  shows  the  modeled  paleogeography  for  the  Channel 
Islands  at  20  ka  and  from  15  ka  to  10  ka  computed  using  either  an 
assumption  of  a  eustatic  sea-level  change  (Fig.  3a)  or  the  gravita¬ 
tionally  self-consistent  ice-age  sea-level  theory  (Fig.  3b).  Fig.  3c 
shows  the  difference  in  the  predicted  shoreline  location  for  each  of 
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Fig.  2.  Regional  maps  of  relative  sea  level  (RSL)  from  20  ka  to  6  ka  predicted  using  the  gravitationally  self-consistent  sea-level  theory  and  the  ANU/LM  model  of  GIA  (see  text). 
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Fig.  3.  The  paleogeography  for  the  Channel  Islands  region,  CA,  at  20  ka  and  from  15  to  10  ka  (at  1-kyr  intervals)  computed:  (a)  with  an  assumption  of  a  eustatic  history:  and  (b) 
when  accounting  for  the  relative  sea  level  (RSL)  history  shown  in  Fig.  2.  (c)  The  difference  in  area  of  exposed  shelf  for  the  Channel  Islands  region  computed  using  the  eustatic  history 
(frame  a)  and  gravitationally  self-consistent  RSL  history  (frame  b).  The  areas  shown  in  blue  represent  zones  exposed  in  the  eustatic  predictions  but  not  in  the  GIA  calculation, 
whereas  areas  in  red  are  zones  exposed  in  the  latter  but  not  the  former.  (For  interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of 
this  article.) 


these  snapshots.  The  predicted  paleoshoreline  position  based  on  a 
eustatic  sea  level  assumption  can  be  in  error  by  several  km  at,  for 
example,  12  ka,  corresponding  to  the  time  of  early  occupation  of  the 
islands  (Erlandson  et  al.,  1996,  2005;  Rick  et  al.,  2001,  2005).  We 

conclude  that  the  sea-level  history  near  the  Channel  Islands  is 
significantly  impacted  due  to  GIA  even  though  the  region  is 
-1500  km  from  the  nearest  North  American  ice  sheets. 

4.  Oregon-Washington  continental  shelf 

Most  of  the  recorded  archaeological  sites  along  coastal  regions 
of  Oregon  and  Washington  are  younger  than  5000  years  ago 
(Lyman,  1991),  with  only  two  older  sites  identified:  the  Indian 
Sands  site  on  the  southern  Oregon  coast  (Davis  et  al.,  2004; 
Erlandson  et  al.,  2008)  and  the  Manis  site  in  northern 


Washington  (Waters  et  al.,  2011b).  Some  have  argued  that  the  lack 
of  older  sites  suggests  late  occupation  of  the  region  (Lyman,  1991; 
Ross,  1990),  but  it  more  likely  reflects  coastal  erosion  or  submer¬ 
gence  of  offshore  sites  associated  with  post-glacial  sea-level  rise 
amplified  by  periodic  subsidence  earthquakes  and  tsunamis 
(Erlandson  et  al.,  2008). 

The  relation  between  local,  site-specific  RSL  change  for  the 
continental  shelf  off  of  southern  Oregon  (42.5°N)  and  the  eustatic 
sea  level  curve  is  similar  to  that  described  above  for  the  Channel 
Islands,  with  differences  of  up  to  -20  m  and  RSL  reaching  present- 
day  levels  -6000  years  later  than  predicted  under  the  eustatic 
assumption  (Fig.  4).  As  one  moves  north  to  the  shelf  off  central 
Oregon  (44.5°N),  however,  RSL  is  below  eustatic  throughout  the  full 
interval,  with  differences  of  as  much  as  -50  m  at  -12  ka  (Fig.  4). 
Further  north  on  the  shelf  off  southern  Washington  (46.5°N),  the 
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Fig.  4.  Examples  of  relative  sea  level  curves  (blue)  for  four  sites  along  the  Ore- 
gon-Washington  coast  compared  to  the  eustatic  curve  associated  with  the  ANU 
deglaciation  history  (red).  (For  interpretation  of  the  references  to  color  in  this  figure 
legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 


relation  between  RSL  and  eustatic  sea  level  is  again  like  that  on  the 
shelf  off  southern  Oregon,  with  RSL  initially  above  eustatic  until 
~15  ka,  when  the  sign  of  the  discrepancy  changes  and  the  predicted 
RSL  is  as  much  as  ~40  m  lower  than  eustatic.  Closer  to  the  margin  of 
the  Cordilleran  Ice  Sheet  in  northern  Washington,  this  relation  then 
changes  significantly  as  the  large  isostatic  and  gravitational  effects 
of  the  nearby  ice  sheet  begin  to  dominate  the  signal.  For  example,  at 
a  site  at  48.5° N,  RSL  was  50-80  m  above  present  level  from  20  to 
15  ka  (Fig.  4),  suggesting  that  sea-level  rise  associated  with  the 
addition  of  meltwater  was  offset  by  sea-level  fall  due  to  isostatic 
uplift  (or  rebound)  of  the  crust  and  the  gravitationally-driven 
migration  of  water  away  from  the  melting  ice  sheet.  From  15  to 
12  ka,  the  latter  mechanisms  continue  to  dominate  the  meltwater 
addition  component  and  sea  level  is  predicted  to  fall  -75  m.  After 
12  ka,  the  eustatic  signal  associated  with  meltwater  addition 
dominates,  and  sea  level  rose  toward  its  present  value.  Over  this 
time  interval,  the  hinge  point  separating  regions  of  post-glacial 
uplift  from  sites  experiencing  peripheral  bulge  subsidence 
migrated  in  this  particular  simulation,  and  therefore  bulge  subsi¬ 
dence  also  contributes  to  the  predicted  sea-level  rise.  Finally, 
further  north  in  coastal  areas  of  Canada  that  were  formerly  covered 


by  the  ice  sheet,  a  monotonic  fall  in  RSL  over  the  last  20,000  years  is 
associated  with  the  combined  effects  of  crustal  uplift  and  the 
gravitationally-driven  migration  of  water  away  from  the  melting 
ice  sheet. 

Fig.  5a  and  b  shows  the  paleogeography  for  the  Ore- 
gon-Washington  continental  shelf  from  20  ka  to  6  ka  computed 
either  assuming  eustatic  sea  level  change  or  accounting  for  ice-age 
adjustments,  respectively.  Given  the  bathymetry  of  this  region, 
there  are  some  notable  differences  between  these  two  paleogeo- 
graphic  reconstructions.  As  an  illustration,  we  focus  on  two  sites  at 
two  times  that  have  implications  for  predictive  modeling  of 
offshore  sites  (Fig.  6). 

The  first  case  is  for  the  central  Oregon  shelf  at  14  ka.  An 
assumption  of  eustatic  sea  level  change  predicts  a  series  of  small 
islands  at  this  time  as  most  of  the  shelf  was  submerged  (Fig.  6a).  In 
contrast,  the  prediction  based  on  the  computed  RSL  history  shows 
that  much  of  the  shelf  in  this  area  remained  emergent,  with  an 
extensive  coastline  (Fig.  6b).  The  second  case  is  for  the  area  close 
to  and  just  north  of  where  the  Columbia  River  enters  the  Pacific. 
An  assumption  of  eustasy  predicts  that  much  of  the  estuarine 
regions  along  this  part  of  the  coast  were  submerged  at  6  ka 
(Fig.  6c),  whereas  the  prediction  using  the  RSL  calculation  suggests 
that  the  shoreline  remained  seaward  of  the  present  coast.  These 
estuaries  are  predicted  to  have  remained  emergent  at  this  time, 
and  likely  would  not  have  served  as  coastal  occupation  sites 
(Fig.  6d). 

Fig.  5c  shows  regional  maps  of  the  predicted  RSL  from  20  ka 
to  6  ka  using  the  gravitationally  self-consistent  sea-level  theory. 
As  with  the  Channel  Islands,  these  maps  show  that  relative  sea- 
level  fall  was  not  spatially  uniform,  but  in  comparison  to  Fig.  2, 
there  are  now  much  stronger  gradients  in  RSL  associated  with 
GIA  due  to  the  greater  proximity  to  the  Cordilleran  Ice  Sheet.  In 
particular,  there  is  a  strong  SSW-NNE  gradient  in  RSL  on  the 
northern  Washington  shelf,  especially  from  20  to  15  ka.  As 
illustrated  by  Fig.  4,  RSL  in  this  northern  region  over  this  time 
interval  was  50-80  m  above  present  at  the  LGM,  indicating  that 
any  sites  from  that  time  period  would  be  exposed  today  along 
the  coast. 

As  a  final  point  we  emphasize  that  the  RSL  predictions  in 
Figs.  4-6  are  based  on  a  single  viscoelastic  Earth  model  (LM)  and 
ice  history.  The  discrepancies  between  these  predictions  and 
eustatic  sea  level  will  be  a  function  of  the  adopted  Earth  model  and 
ice  history,  but  predictions  based  on  all  plausible  Earth  models  will 
show  comparable  levels  of  discrepancy. 

5.  Bering  Sea  continental  shelf 

During  the  LGM  sea-level  lowstand,  nearly  2  x  106  km2  of 
shallow  continental  shelf  underlying  the  Bering  and  Chukchi  Seas 
was  exposed,  creating  the  contiguous  land  mass  referred  to  as 
Beringia  that  extended  from  eastern  Siberia  east  to  the  Mackenzie 
River.  The  significance  of  this  land  bridge  as  a  migration  route  be¬ 
tween  Asia  and  North  America  for  animals  and  early  people  has 
long  been  recognized  (Hopkins,  1967, 1973),  and  constraining  when 
it  became  submerged  is  important  to  understanding  when  such 
land  migrations  ended.  The  timing  is  also  important  for  assessing 
the  influence  of  the  opening  of  a  water  connection  between  the 
Pacific  and  Arctic  Oceans  on  regional  climate  and  ocean  circulation 
(Hu  et  al.,  2010). 

Recognizing  that  the  concept  of  eustatic  sea  level  change  should 
not  be  applied  to  reconstruct  Beringian  paleogeography,  McManus 
and  Creager  (1984)  constrained  regional  paleo  water  depths  using 
marine  faunal  estimates  and  concluded  that  the  initial  submer¬ 
gence  of  the  land  bridge  occurred  -17,500  years  ago.  In  contrast, 
Elias  et  al.  (1996)  dated  terrestrial  plant  remains  from  cores  taken 
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Fig.  5.  The  paleogeography  and  relative  sea  level  (RSL)  history  for  the  Oregon-Washington  continental  shelves  from  20  ka  to  6  ka.  (a)  The  paleogeography  predicted  using  an 
assumption  of  a  eustatic  sea-level  change,  (b)  The  paleogeography  predicted  by  accounting  for  the  RSL  history  in  Fig.  5c.  (c)  Regional  maps  of  RSL  predicted  using  the  gravitationally 
self-consistent  sea-level  theory  and  the  ANU/LM  GIA  model. 


on  the  Chukchi  and  Bering  Sea  shelves  that  suggested  that  the 
submergence  did  not  occur  until  -13,000  years  ago,  and  argued  that 
the  older  ages  from  McManus  and  Creager  (1984)  likely  were 
contaminated  by  old  carbon. 

Fig.  7  shows  our  sea  level  predictions  at  each  of  the  three  sites 
where  Elias  et  al.  (1996)  obtained  constraints  on  RSL.  Each  frame 
includes  the  eustatic  curve  for  the  ANU  ice  history  (red  line)  and 
GIA  predictions  based  on  the  ANU/LM  ice  and  Earth  model  pairing 
(blue  line).  Clearly,  an  accurate  prediction  of  the  space  and  time 
dependence  of  the  inundation  of  Beringia  requires  that  one  account 
for  ice-age  adjustments  on  sea  level.  Moreover,  RSL  predictions 
based  on  the  ANU/LM  GIA  model  in  Fig.  7  are  consistent  with  the 
RSL  constraints  provided  by  the  limiting  ages  on  terrestrial  organics 


from  Elias  et  al.  (1996),  insofar  as  the  presence  of  the  organics 
suggests  that  they  were  not  submerged  until  sometime  after  they 
lived. 

Fig.  8a  shows  the  paleogeography  of  the  Bering  Sea  conti¬ 
nental  shelf  from  15  ka  to  6  ka  under  the  assumption  of  a  eustatic 
sea-level  history  (note  that  there  are  no  substantial  differences  in 
the  paleogeography  of  this  region  between  20  ka  and  15  ka). 
Fig.  8b  shows  the  paleogeography  for  the  Bering  Sea  continental 
shelf  from  15  ka  to  6  ka  computed  using  the  RSL  history  in  Fig.  8c. 
There  are  some  notable  differences  between  the  two  re¬ 
constructions.  In  particular,  we  note  three  examples  that  have 
implications  for  predictive  modeling  of  offshore  sites  as  well  as 
for  the  use  of  the  Bering  land  bridge  for  migration.  First,  both 
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reconstructions  show  an  elaborate  archipelago  to  the  south  of 
the  Bering  land  bridge  that  may  have  served  as  occupation  sites, 
but  with  significant  differences  at  the  local  scale  that  will  be 
important  to  account  for  in  predictive  modeling.  Second,  the 
assumption  of  eustasy  predicts  that  the  Norton  Sound  (at  bottom 
right  on  each  frame)  remained  largely  emergent  until  -12  ka;  in 
contrast,  inundation  of  the  Sound  is  predicted  to  have  begun 
-1000  years  earlier  at  -13  ka  when  GIA  effects  are  incorporated 
into  the  reconstruction.  Third,  in  the  paleogeographic 


reconstruction  based  on  eustasy,  inundation  of  the  northern 
Bering  land  bridge  in  the  Chukchi  Sea  region  began  at  -11  ka, 
whereas  the  RSL  calculation  predicts  inundation  -10  ka, 
demonstrating  that  GIA  effects  had  a  significant  impact  on  the 
evolution  and  ultimate  disappearance  of  the  land  bridge.  Lastly, 
Fig.  8c  shows  the  predicted  RSL  history  from  15  ka  to  6  ka  based 
on  the  ANU/LM  GIA  model.  The  regional  predictions  exhibit  a 
spatial  gradient  of  -30  m  from  southwest  to  northeast  at  15  ka 
and  this  gradient  largely  disappears  by  -10  ka. 
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Fig.  6.  Paleogeography  of  the  central  Oregon  continental  shelf  at  14  ka  predicted  assuming  either  (a)  a  eustatic  history,  or  (b)  a  gravitationally  self-consistent  RSL  history  based  on 
the  ANU/LM  GIA  model,  (c,  d)  Analogous  to  frames  (a,  b)  except  for  the  continental  shelf  near  the  mouth  of  the  Columbia  River  at  6  ka. 
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Fig.  7.  Predictions  of  relative  sea  level  change  for  three  sites  within  the  Beringian 
region  based  on  the  ANU/LM  (blue  line)  GIA  models  (see  text).  Also  shown  is  the 
eustatic  curve  associated  with  the  ANU  deglaciation  history  (red  line).  The  diamonds 
on  each  frame  represent  RSL  constraints  (see  text)  recovered  from  marine  cores  (Elias 
et  al.,  1996).  The  three  sites  are  located  at:  (a)  65.17°N,  167.5°W,  (b)  64.08°N,  164.5°W, 
and  (c)  71°N,  166°W.  (For  interpretation  of  the  references  to  color  in  this  figure  legend, 
the  reader  is  referred  to  the  web  version  of  this  article.) 


Clark  et  al.  (2014)  considered  predictions  of  GIA-induced 
Beringian  sea  level  change  using  a  different  ice  history/Earth 
model  pairing;  specifically,  the  VM2  radial  profile  of  mantle  vis¬ 
cosity,  which  is  coupled  to  the  ICE-5G  ice  history  since  the  last 
interglacial  (Peltier,  2004).  The  VM2  model  is  characterized  by  a  90- 
km  thick  elastic  lithosphere  and  a  moderate,  factor  of  ~4-5  increase 
in  viscosity  from  the  upper  mantle  (~5  x  1020  Pa  s)  to  the  lower 
mantle  (~2  x  1021  Pa  s). 

Fig.  9a  shows  the  paleogeography  of  the  Bering  Sea  continental 
shelf  from  15  ka  to  6  ka  under  the  assumption  of  a  eustatic  sea-level 
history  based  on  ICE5G.  Fig.  9b  shows  the  paleogeography  for  the 
Bering  Sea  continental  shelf  from  15  ka  to  6  ka  based  on  the  RSL 
history  predicted  by  adopting  the  ICE5G/VM2  GIA  model  (Fig.  9c). 
We  focus  on  comparing  the  paleogeographies  computed  using  the 
two  GIA  predictions  (Fig.  8b  compared  to  Fig.  9b).  There  are  two 
notable  differences  that  have  relevance  to  the  history  of  the  Bering 
land  bridge  and  possible  coastal  migration  routes  and  occupation 
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Fig.  8.  Predicted  paleogeography  and  relative  sea  level  (RSL)  history  for  the  Beringian  continental  shelf  from  15  ka  to  6  ka.  (a)  Paleogeography  computed  using  an  assumption  of  a 
eustatic  sea  level  change  derived  from  the  ANU  deglaciation  history  (Fleming  et  al.,  1998).  (b)  Paleogeography  when  accounting  for  the  gravitationally  self-consistent  RSL  changes  in 
(c),  where  the  latter  was  computed  using  the  ANU/LM  GIA  model. 


sites.  First,  the  RSL  prediction  based  on  the  ICE5G/VM2  GIA  model 
shows  that  initial  submergence  of  the  Bering  Straits  occurred 
-11  ka,  or  -1000  years  earlier  than  the  prediction  based  on  the 
ANU/LM  GIA  model.  Second,  the  ICE5G/VM2  GIA  model  results 
predict  that  much  of  the  Norton  Sound  remained  submerged  since 
the  LGM,  in  contrast  to  the  ANU/LM  GIA  model,  which  predicts  that 
the  inundation  of  the. 

6.  A  sensitivity  analysis 

The  predictions  described  thus  far  were  based  on  1-D  viscoelastic 
Earth  models;  specifically,  model  LM  in  Figs.  1-8  and  VM2  in  Fig.  9. 
However,  the  west  coast  of  North  America  is  characterized  by  an 


active  plate  boundary  and  variations  in  both  lithospheric  thickness 
(Watts,  2001 )  and  mantle  structure  (Ritsema  et  al.,  2011 ).  To  consider 
the  impact  that  this  complexity  might  have  on  the  predictions  of  sea- 
level  change  and  shoreline  migration,  we  ran  a  sea-level  simulation 
on  a  finite  volume  numerical  model  of  the  GIA  process  that  is  capable 
of  incorporating  such  features  (Latychev  et  al.,  2005).  The  3-D  Earth 
structure  adopted  in  the  simulation  is  described  in  detail  in 
Austermann  et  al.  (2013 ),  although  in  the  present  case  we  incorporate 
lateral  variations  in  viscosity  throughout  the  mantle  and  superim¬ 
pose  these  on  the  1-D  LM  viscosity  profile.  The  variations  reach  five 
orders  of  magnitude  peak-to-peak. 

Fig.  10  shows  predictions  of  relative  sea  level  change  at  the  four 
sites  along  the  Oregon-Washington  coast  considered  in  Fig.  4.  The 
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Fig.  9.  Predicted  paleogeography  and  relative  sea  level  (RSL)  history  for  the  Beringian  continental  shelf  from  15  ka  to  6  ka.  (a)  Paleogeography  computed  using  an  assumption  of  a 
eustatic  sea  level  change  derived  from  the  ICE5G  deglaciation  history  (Peltier,  2004).  (b)  Paleogeography  when  accounting  for  the  gravitationally  self-consistent  RSL  changes  in  (c), 
where  the  latter  was  computed  using  the  ICE5G/VM2  GIA  model  (Peltier,  2004). 


blue  lines  are  predictions  based  on  the  1-D  viscoelastic  Earth  model 
ANU  (i.e.,  these  results  are  reproduced  from  Fig.  4),  while  the  red 
lines  are  computed  using  the  3-D  viscoelastic  Earth  model.  Both 
simulations  adopt  the  ANU  ice  history.  The  perturbations  to  the 
predictions  associated  with  lateral  variations  in  structure  are  non¬ 
monotonic  in  time  and  in  very  general  terms  the  3-D  predictions 
are  shifted  earlier  by  -1  kyr  relative  to  the  1-D  simulation.  (The 
result  for  the  site  at  46.5°N  at  times  prior  to  15  ka  is  a  clear 
exception.)  We  conclude  that  in  archaeological  applications  where 
predictions  of  shoreline  evolution  require  a  precision  better 
than  ~  1  kyr,  the  sea-level  simulations  adopted  in  such  analyses 
should  incorporate  the  full  complexity  of  Earth  structure  beneath 
western  North  America. 


7.  Conclusions 

The  evidence  for  pre-Clovis  occupation  of  the  Americas  as  early 
as  -15-16  ka  provides  additional  support  for  the  Fladmark  (1979) 
hypothesis  that  the  first  Americans  migrated  from  Asia  to  lands 
south  of  the  North  American  ice  sheets  by  a  route  along  the  west 
coast  of  North  America.  Constraining  possible  migration  pathways, 
as  well  as  identifying  potential  early  occupation  sites  that  are  now 
below  sea  level,  requires  an  accurate  retrodiction  of  sea-level 
change  during  the  last  deglaciation  and  through  the  Holocene. 
We  have  shown  that  departures  from  a  uniform  (eustatic)  deglacial 
sea-level  rise  due  to  deformational,  gravitational,  and  rotational 
effects  driven  by  the  exchange  of  mass  between  ice  sheets  and 
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Fig.  10.  Predicted  relative  sea  level  curves  for  four  sites  along  the  Oregon-Washington  coast.  The  blue  lines  are  predictions  based  on  the  1-D  viscoelastic  Earth  model  LM  and  the 
ANU  deglaciation  history  (these  lines  are  reproduced  from  Fig.  4).  The  red  lines  are  analogous  to  the  blue  lines  except  that  the  predictions  are  based  on  the  3-D  viscoelastic  Earth 
model  described  in  the  text.  (For  interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 


oceans  during  the  last  deglaciation  (i.e.,  glacial  isostatic  adjust¬ 
ment)  led  to  substantial  impacts  on  the  paleogeography  across  the 
California-Oregon-Washington  and  Bering  Sea  continental 
shelves.  These  departures  must  be  accounted  for  in  any  paleogeo- 
graphic  reconstruction  of  these  regions.  Numerical  predictions  of 
post-glacial  sea-level  changes  described  herein  were  based  on  a 
specific  pairing  of  ice  history  and  Earth  model.  The  paleogeographic 
reconstructions,  however,  will  be  sensitive  to  uncertainties  in 
either  input.  In  addition,  our  modeling  did  not  include  the  impact 
on  sea  level  of  sediment  redistribution,  which  may  be  significant  in 
the  vicinity  of  the  Columbia,  Eel  and  Klamath/Trinity  river  systems. 
Further  work  is  required  to  address  each  of  these  issues  before 
definitive  predictions  of  the  location  of  now-submerged  archaeo¬ 
logical  sites  are  possible. 
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